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Motivated by excellent comparisons of spin wave theory (SWT) with Quantum Monte Carlo 
(QMC) calculations of non-interacting hardore bosons (Bernardet et al.) and superlattices (Hen et 
al.), we develop an SWT for interacting hardcore bosons on a d-dimensional cubic lattice at zero 
temperature and compare thermodynamical properties with results of Stochastic Series Expansions 
(SSE): there is good comparison between SWT and SSE results for the simple cubic lattice at small 
repulsive strengths but differences arise between the results of the two approaches for the square 
lattice at strong repulsions. The boundaries of the three phases - Mott insulating (MI), superfluid 
(SF) and Neel solid (NS) - can be readily estimated from calculations of particle density, signalling 
either a first order (SF -H- NS) or second order transition (MI -H- SF). Comparison is made with 
results of the Bethe ansatz in one dimension and QMC in two dimensions. 



Bosons on a lattice appear as a natural model for 
several physical systems. Such a model was used to elu- 
cidate the A-transition in liquid helium Ili which, at the 
time, was still lacking a microscopic basis. Furthermore, 
the authors in Ref. [1] showed that the partition function 
of a suitably approximated model - which later came to 
be variously described as the hardcore boson Hubbard 
model (HCB) - is equal to the partition function of a 
system of vector spins - the Hciscnbcrg model in a longi- 
tudinal magnetic field - as defined on an infinite lattice, 
when a suitable identification of constants and operators 
between the two models is made. The latter system's 
partition function, they argued, could be calculated by 
an approximation scheme like spin- wave analysis. More- 
over, an additional reason for preferably treating the cor- 
responding spin model is one of calculational convenience 
vis-a-vis keeping track of the unusual commutation rela- 
tion of hardcore bosons (see next section). 
Within this approximation, several studies at zero and fi- 
nite temperatures^— were undertaken. Although Ref. [2] 
focussed on describing the phase boundaries, many ther- 
modynamical quantities were not explicitly calculated. 
The authors of Ref. [3] repeated the calculations, on a 
square lattice mostly near half-filling, in order to resolve 
some previously conflicting results regarding the spec- 
trum in the earlier predicted phases. The focus of the 
present work, unlike previous reports, will not be on the 
existence or otherwise of a supersolid phase but rather on 
calculating thermodynamic quantities using a collinear 
spin wave theory (SWT) in the interacting model and 
comparing the results with other methods. 
In fact, with the subsequent development of sophisticated 
Quantum Monte Carlo (QMC) techniques, SW treatment 
has been seen to provide a surprisingly adequate descrip- 
tion for many physical observables in the HCB model like 
particle density, ground state energy and suchlike for (a) 
the non-interacting case^ and (b) superlattice systems^. 
Indeed, as Bernardet et al. and Hen et al. showed, 
the discrepancies, if any, between the results of the two 



approaches - QMC and SWT - are extremely small for 
many thermodynamical properties. This provides further 
incentive to revisit the corresponding nearest neighbour 
(NN) interacting problem for a general cubic lattice us- 
ing SWT. 

The goal of the paper is to provide an SW treatment for 
the interacting HCB model on a d-dimensional cubic lat- 
tice, calculate certain thermodynamic quantities of inter- 
est and compare our results with some known solutions, 
where available, in any particular dimension. In fact, 
it was suggested^ that SWT should be sufficiently good 
even in the interacting case; we find that, starting from 
a collinear phase, the comparison between SWT and re- 
sults of Stochastic Series Expansion (SSE) is reasonable 
but not as good as for the non-interacting case for the 
square lattice with strong repulsions V. For the simple 
cubic lattice, there is good agreement between the two 
approaches at a lower V value. 



I. MODEL 

The hardcore boson model with nearest neighbour 
interaction on a d-dimensional cubic lattice is represented 
by the Hamiltonian 



<ij> 



<ij> 



The 



annihilation and creation field operators 



satisfy the usual bosonic commutation relations for 
different sites and fermionic anticommutation relations 
for the same site, thereby preventing double occupancy; 
h = a'a is the boson number operator. The interaction 
V and hopping t are between NN sites, and the chemical 
potential fj, controls the particle density. The spinlcss 
bosons can be mapped to an anisotropic Hciscnbcrg 
spin-1/2 model with the following equivalence relations: 
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site i. 



where the mean field particle density 



II. SPIN WAVE ANALYSIS 

The mean field and spin wave analyses proceed in the 
familiar textbook way^. In the subsequent formulation, 
we closely follow the notation of Rcf. [5] . 
Upon changing to the spin representation, the Hamilto- 
nian in Eq. (1) becomes 



H 



<ij> 



<ij> 



-0i-zV/4)N, (2) 



where z is the number of nearest neighbours and N 
is the number of lattice points. For cubic lattices in 
e?-dimcnsions, z = 2d. 

Consider the mean- field ground state \tp) = 
niLi(sin(0/2) + cos(0/2)a|)|O), where 9 is the an- 
gle made by the (pseudo)spins with respect to the z-axis. 
We choose, without loss of generality, an ordering such 
that the xz plane contains all the spins. Minimizing the 
grand canonical potential per site f2 ~ {tp\'H\ip) with 
respect to the angle 9, we obtain, for the mean field 
values, 



^MF = PMF{-^VpMF - P + ztpMF 



zt), 



Pmf = cos 2 (9/2) 



zt 



z(2t + vy 



The energy density is given by Emf = ^mf + PPmf- 
Spin wave corrections can be included by twisting the 
2-axis to align with the direction of MF magnetization: 

Sf = Sf cos (9) + §? sin (6) 
Sf = -§f sin (6) + Sf cos (9). 

(3) 

In this rotated frame, the excitations are bosonic in na- 
ture because the spin-flips change the S z projected spin 
values by ±1. This requires a bosonic representation us- 
ing the Holstcin-Primakoff transformations for hardcore 
bosons 



S : 



Sf 



1 



i 



A) 
-h) 



(4) 



With this substitution the spin wave Hamiltonian, up to 
quadratic terms, is readily shown to be 



E {( 6 H + W(| + T) sin2 V) + ( b h + + cos 2 (9)) + \ sin 2 (9))} 
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E b l b i t^ sin2 W ~ T cos2 W + - ^) cos (9)} +N[(^- |) cos (9) + cos 2 (9) t sin 2 (9)) 
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E[ M f /2 sin w -j( 2t + v ) cos ( ) shi (*)] $ + b <)- 



(5) 



So that the particle number is conserved in the new 
bosonic system, we require that each of the last sum- 
mands disappear. This gives the spin-wave condition for 
9 to be 



cos (9) 



2/i - zV 
z(2t + V)' 



This was the same condition previously obtained with- 
out the spin waves as well, noted as well in the non- 
interacting case£. We mention in passing that we 
can, from the above, obtain a d-dimcnsional analog of 
the Pokrovsky-Talapov linei^ separating the massless 
and the ferromagnetic phases in the one-dimensional 



anisotropic Heisenberg model in a longitudinal field by 
equating cos (9) = 1. This gives the second order phase 
separating line between the MI and SF phases to be 



{p) p c T = z(t + V). 



(6) 



The spin wave Hamiltonian retains the same form as 
in Ref. [5] except for the coefficients of the operators. 
Therefore the Bogoliubov transformation 



be = ucO-c — v.-a - 

k k k k _k 



bl = u - u k«_k- 
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k k 
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FIG. 1: (Colour online) (a) Mean field particle density on the square lattice compared with spin-wave results and implementation 
of Stochastic Series Expansion using the ALPS library^ at V/t = 10. Close to the first-order transition to the checker-board 
solid, the SWT does not converge. Inset shows similar comparisons of boson density for the simple cubic lattice at V/t = 1. (b) 
Energy density for the simple cubic lattice in the full boson density range, with arrows indicating MF values at p = 1. Inset 
shows the boson density computed using spin-wave and ALPS's SSE* on the square lattice at V = 0, verifying the result of 
Ref. [5]. 



to separate the MF ground state and the excited states 
above it - represented by the a bosons - immediately 
simplifies to Eq. (22) in Ref. [5] but with different ex- 
pressions for the coefficients i.e. 



W-sw = H 



MF 



E 

k 



The sum over k is over the first Brillouin zone and the 
mean field Hamiltonian corresponds to the mean field 
grand canonical potential £Imf- The variables A^B^ 
are here given by 



- ' ^[(V-2t) 



Ac = 
k 2 4 



{2n-zVf 
z 2 {2t + V) 



B~ = ^k 

C z 2 (2t + V) 



[ti?(t + V) + n(zV-n)], 



7k = E cos 

i=l 

Here the lattice constant was taken to be one and gauge 
invariance in the a operators is ensured by choosing 



1 



A * ±1). 



III. THERMODYNAMIC QUANTITIES 

Using the above results, we can generalize the results 
in Ref. [5] to a cubic lattice in arbitrary dimensions. 



The average particle density, including the spin waves 
but calculated within the MF state, is given by 



Psw 



sir 



2/x - 



Pmf 



zV 



'z 2 N(2t + V) ^ 

k 



> £ ( Jk=£g= _ 1), (9) 



Al 



Bl 

k 



where £lsw is given by the expectation value of the first 
line of Eq. (8), whereas the coefficient — B\ m tnc 

second line represents the dispersion of the excitations. 
In Fig. 1(a) we compare the boson density calculated 
using SWT, with results of Quantum Monte Carlo (SSE 
implemented using the ALPS library^) simulations and 
MF analysis on the square lattice for V/t = 10. In the 
SSE simulation we used /3 = 1/T = 40 on a 25x25 system 
using half a million sweeps and a thermalization of ten 
thousand; the error bars are not shown (because they are 
only as large as the symbol size) and simulations were 
done for various system sizes (L = 25,12,6,4), sweeps, 
thermalization and inverse-temperatures ((3 = 40, 20, 10) 
to check for convergence. As can be seen SWT is 
an improvement over the MF results but compares 
only approximately with SSE, with the differences 
increasing towards the phase boundary; the comparison 
is not as good as with the non-interacting system in 
Ref. [5], which we have verified as shown in the inset 
of Fig. 1(b). For one dimension higher i.e. simple 
cubic lattice, and smaller V/t, spin wave analysis 
compares well with an SSE simulation using the ALPS 
library as shown in the inset of Fig. 1(a) for V/t = 1; 
the system size used was 5x5x5 with the same ther- 
malization, sweeps and temperature as mentioned above. 
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FIG. 2: (Colour online) (a) Quasiparticle - a bosons - dispersion for 2 interaction strengths along k x — k y = k z in the 3D 
simple cubic lattice. As explained in the text, ferromagnetic and antiferromagnetic ordering is readily seen from behaviour 
near the zone edge or centre, (b) SWT phase diagram compared with other known results like QMC— on a square lattice and 
Bethe ansatz on a linear chain. Naturally, we expect the SWT to be more accurate in higher dimensions. 



The energy density in the SW approximation is ob- 
tained as before within mean field i.e. E$w = ^sw + 
ppsw- Finally, the functional expressions for the ex- 
pectation value of the condensate fraction p - particle 
density in the zero momentum sector - within MF and 
SW remain the same as in Rcf. [5] in the absence of NN 
interactions: 

P0,MF = PMfQ — PMF), 

Po,sw = p(l - p) - j^^2 v l- (10) 

The energy density in the simple cubic lattice for 2 differ- 
ent interaction strengths V are computed as functions of 
the SW particle density by Legendrc transforming from 
the p variable and shown in Fig. 1(b). The condition 
E(p) = E(l — p) is not satisfied at a finite V because it is 
clear that larger particle densities contribute more to the 
intersite repulsive energy. From the figure, the usual ex- 
pectations are borne out: higher intersite repulsion adds 
more SW contributions to the MF value. That is, lower 
the interaction, better is the MF description. This may 
be seen by the distance between the location of the ar- 
rows (MF values) and the corresponding SW values for 
two different interaction strengths. It is worth noting 
from Eq. (10) that the MF condensate fraction depends 
only on the particle density and not on the interaction 
strength; whereas the larger the V value, the more will 
Bose-Einstein condensation be suppressed within SWT. 



IV. EXCITATIONS 

As seen in the previous section, the excitation spec- 
trum of the a bosons is given by U)(k) = ./A! — B"l. 



The quantities Ag, -Bg still depend on p,V but the de- 
pendence is notationally suppressed. The dispersion is 
calculated along the line k x = k y = k z for a simple cubic 
lattice and plotted in Fig. 2(a) for two different interac- 
tion strengths. 

Far away from the expected checkerboard solid phase, 
antiferromagnetic (AF) ordering is not expected for any 
particle density. This is indeed true and checked for 
V/t = 1 for all range of densities between the empty and 
full phases. That is, the dispersion in this case is always 
gapped at k = (tt, 7T,7r). The dispersion at the bound- 
aries of these two ferromagnetic phases - p/V = —6, 12 
- show limigi^g w(k) oc |k| 2 (upper panel of Fig. 2(a)) 
as expected for such ordering. The gap to magnon ex- 
citations in the MI phase is of course just given by 
A = (M-(M)f T )/2. 

In the strong coupling limit - V/t — 10 - antiferromag- 
netic ordering sets in as the checkerboard solid is ap- 
proached. This is seen in the lower panel of Fig. 2(a), 
where for p/V = 0.6, 

\|k- (tt, 7T, 7r)|, |k| « (7r,7r,7r). 

Finally, although the SW ground state phase dia- 
gram has already been fleshed out in detail^, our pur- 
pose here is merely to compare it with other known re- 
sults. The SW construction proceeds in the usual way: 
for every interaction strength, the particle density is cal- 
culated for all chemical potentials, p will have a discon- 
tinuity at two values of p corresponding to the first order 
phase transition between the solid and superfluid. These 
will define the boundaries of the expected checkerboard 
solid lobe. Moreover, it must be recalled that within our 
SWT, only collincar phases are describable because we 
chose 6i — 6. Therefore, we do not expect convergence 
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in the solid phase anyway. With these assumptions, the 
SWT phase diagram is constructed for the HCB model on 
any d-dimensional cubic lattice; the separation between 
the massless and the ferromagnetic phase was already 
given by Eq. (6). 

This is compared to QMC calculations in two 
dimensions^ and the Bethe ansatz solution in one dimen- 
sion. The latter was solved for the anisotropic Hciscnberg 
model on a chain-^; the critical line (y) c separating the 
massless and antiferromagnetic phase is given by 

-[(-^) c -l]=sinh(A) ]T (-lrcosh-^nA), (11) 

n— — oo 

where A = cosh -1 (^). The Pokrovsky-Talapov line re- 
mains the same. The comparisons between Bethe ansatz 
in ID, QMC in 2D and SWT in d-D is shown in Fig. 2(b). 

To summarize, we have generalized the SWT for cu- 
bic lattices in arbitrary dimensions for a hardcore boson 
model with nearest neighbour interaction. Based on ear- 
lier work££, where in Ref. [5] it was suggested that such 
an SWT should give reasonable results for the interacting 



system as well, we find that the comparison with SSE is 
not as good as it is with the non-interacting case for the 
square lattice and strongly repelling nearest neighbour 
bosons. Although we did not start with a bond-ordered 
phase, as suggested in Ref. [5], we can compute observ- 
ables like energy density and condensate fraction as long 
as the system is not within the checkerboard solid phase: 
the collinear MF states should, in principle, describe the 
SF and MI phases. These quantities are evaluated for 
the simple cubic lattice and compared with the mean 
field values. In a simple cubic lattice and away from 
the strongly-repulsive regime, we find good comparison 
for the particle-density between results of spin-wave and 
SSE simulations. 
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